{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "96bdf66c-eb7f-40e1-a309-2651ff9dda38",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "import scipy.integrate as integrate\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "f282c1bc-05a3-4674-a236-2f0f1593a92b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Brec_fn(Btype,Bamb,x,PSE):\n",
    "    '''\n",
    "    type: \n",
    "        0 = Constant field, non-spinning chondrule\n",
    "        1 = Constant field, spinning condrule\n",
    "        2 = Sinusoidally varying field, non-spinning chondrule\n",
    "        3 = Sinusoidally varying field, spinning chondrule\n",
    "    '''\n",
    "    if Btype == 0:\n",
    "        #here PSE is the paleointensity standard error\n",
    "        return np.exp(-(x-Bamb)**2 /(2*(PSE)**2))/(np.sqrt(2*np.pi)*PSE)\n",
    "    \n",
    "    elif Btype == 1:\n",
    "        return 1/Bamb*np.ones(len(x))\n",
    "    \n",
    "    elif Btype == 2:\n",
    "        # not defined at 0 so get as close as possible without introducing other artifacts\n",
    "        #n = np.minimum((x/Bamb)**2,np.ones(x.shape)*(1-1e-4))\n",
    "        return 2/(np.pi*Bamb) * 1/np.sqrt(1-(x/Bamb)**2)\n",
    "    \n",
    "    elif Btype == 3:\n",
    "        # arctanh not defined at 1 so get as close as possible \n",
    "        n = np.maximum((x/Bamb)**2,np.ones(x.shape)*1e-9)\n",
    "        return 2/(np.pi*Bamb) * np.arctanh(np.sqrt(1-n))\n",
    "        \n",
    "    else:\n",
    "        raise NameError(\"Type not recognized\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "46c77a07-c13f-4d76-a1c9-601c669bc687",
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood_table(chond_pints, chond_sigmas, Btype=0, \n",
    "                     Bamb_min=10, Bamb_max=250, Bamb_int=2, dx=0.1, PSE=False):\n",
    "    '''\n",
    "    Produces a table of likelihoods for a chosen Brec function and the given data.\n",
    "    \n",
    "    chond_pints: (N,) np array of chondrule paleointensities \n",
    "    chond_sigmas: (N,) np array of 1-sigma uncertainties in chondrule paleointensities\n",
    "    Brec: function describing recorded field given as type value for Brec()\n",
    "       [0] = Constant field, non-spinning chondrule\n",
    "        1 = Constant field, spinning condrule\n",
    "        2 = Sinusoidally varying field, non-spinning chondrule\n",
    "        3 = Sinusoidally varying field, spinning chondrule\n",
    "    Bamb_min: Minimum ambient field to consider in microtesla [10]\n",
    "    Bamb_max: Maximum ambient field to consider in microtesla [250]\n",
    "    Bamb_int: Ambient field step size to use in calculations in microtesla [2]\n",
    "    PSE: Paleointensity standard error.  This is used for the constant ambient field, no-spin case to approximate the uncertainty on the true value. By default, the function automatically computes the standard error of the dataset. One can manually enter a number, for example a small one, to approximate a Dirac delta function (e.g., true single ambient field).\n",
    "    '''\n",
    "    \n",
    "    N = len(chond_pints)\n",
    "    if not PSE:\n",
    "        PSE = np.std(chond_pints)/np.sqrt(len(chond_pints))\n",
    "    \n",
    "    # loop over Bamb\n",
    "    Bamb_list = np.arange(Bamb_min,Bamb_max+Bamb_int,Bamb_int)\n",
    "    \n",
    "    likelist = np.zeros((N,len(Bamb_list)))\n",
    "    \n",
    "    for i,Bamb in enumerate(Bamb_list):\n",
    "        x = np.arange(0,Bamb+dx,dx)\n",
    "        \n",
    "        # extends range of integration to cover full Bamb range\n",
    "        if Btype == 0:\n",
    "            x = np.arange(0,Bamb+5*PSE+dx,dx)\n",
    "        \n",
    "        # deal with asymptote in sinusoidally varying, non-spinning case\n",
    "        if Btype == 2:\n",
    "            x = np.concatenate([x[:-1],np.linspace(x[-2],x[-1],100)[:-2]])\n",
    "        \n",
    "        # Brec\n",
    "        Brec_x = np.tile(Brec_fn(Btype,Bamb,x,PSE),(N,1)).T\n",
    "        \n",
    "        # PDF\n",
    "        x_t = np.tile(x,(N,1)).T\n",
    "        \n",
    "        PDF_x = norm.pdf(x_t, np.tile(chond_pints,(len(x),1)), \n",
    "                         np.tile(chond_sigmas,(len(x),1)))\n",
    "        \n",
    "        likelist[:,i] = integrate.cumtrapz(Brec_x*PDF_x,x_t,axis=0,initial=0)[-1]\n",
    "    \n",
    "    # assemble output\n",
    "    likelihood = {}\n",
    "    likelihood['table'] = likelist.T\n",
    "    likelihood['total'] = np.nansum(np.nanprod(likelist.T,axis=1))*Bamb_int\n",
    "    likelihood['normalized'] =np.nanprod(likelist.T,axis=1)*Bamb_int/likelihood['total']\n",
    "    likelihood['Bamb'] = Bamb_list\n",
    "      \n",
    "    return likelihood, Brec_x, PDF_x\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34fa955d",
   "metadata": {},
   "source": [
    "## Example 1: Find posterior likelihood function and total likelihood for single model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "ce9af6b7-b32c-40ed-9a0c-2a5c55bc5e73",
   "metadata": {},
   "outputs": [],
   "source": [
    "    '''\n",
    "    Example run with Semarkona (Fu et al. 2014) data. See manuscript text for additional explanation.\n",
    "    '''\n",
    "\n",
    "chond_pints = np.array([37.0, 40.0, 24.5, 11.5, 21.0, 0.0, 0.0])\n",
    "\n",
    "chond_sigmas = np.array([11.0, 13.5, 8.7, 4.7, 5.5, 8.7, 8.7])\n",
    "\n",
    "# Btype = 3 corresponds to sinusoidally varying field, spinning chondrule\n",
    "# It is important to choose Bamb_min and Bamb_max so that the >>0 portions of the distribution are fully covered\n",
    "# In other words, the posterior PDF at Bamb_min and Bamb_max should ~0; otherwise the 'total' variable will not be representative\n",
    "# Choose the Bamb_int variable as a balance between resolution and computation time\n",
    "\n",
    "\n",
    "test1 ,Brec,PDF= likelihood_table(chond_pints, chond_sigmas, Btype=3, \n",
    "                                  Bamb_min=10., Bamb_max=250., Bamb_int=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "efaf20f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.702289177451611e-11"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Explanation of output variables:\n",
    "\n",
    "# test1 contains normalized posterior distributions of paleointensity Bamb:\n",
    "#test1['normalized']\n",
    "\n",
    "# test1 constains an array of Bamb values computed:\n",
    "\n",
    "#test1['Bamb']\n",
    "# test1 contains the total (marginal) likelihood of the utilized model given the dataset:\n",
    "test1['total']\n",
    "\n",
    "#Brec and PDF are used internally by the likelihood_table function and should not be commonly needed. \n",
    "#Brec contains arrays of Brec (recorded field) values for the set of specified Bamb (true ambient field)\n",
    "#PDF contains arrays of normalized probabilities for each chondrule data point at a given Brec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "d4cb1c43-c302-4fd9-91cc-e20a31569429",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fa098c7e2b0>"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWGklEQVR4nO3db4xc1X3G8edha5JNk2aT4lZmTWqnslBpUW20AiRXeYGaYpO2dnlRSNQmRZEs1LgKKKVamhelfYNblJAgISKSoIa2CiCFulagcqI6VVRUEq9j88dxnWwIEV5bwVEwSWQLbPPri7ljxuM7M3dn7sz99/1Ilr1z77Dn7CzP3Dnnd851RAgAUF8XFd0AAMB4EfQAUHMEPQDUHEEPADVH0ANAzf1S0Q1Ic8kll8SaNWuKbgYAVMa+fft+EhEr046VMujXrFmjhYWFopsBAJVh+0e9jjF0AwA1R9ADQM0R9ABQcwQ9ANQcQQ8ANVfKqpsm2rl/SffsPqyjJ07p0plp3XH95dq6YbboZgGoAZdx98q5ubloUnnlzv1LuvPx53Tq9Nlzj1lSSJol9AFkYHtfRMylHWPopgTu2X34vJCXWiEvSUsnTunOx5/Tzv1Lk28YgFog6Au0c/+SNu7Yo6UTp/qed+r0Wd2z+/CEWgWgbgj6grSHawaFfNvSiVNc1QMYCkFfkLThmkEYwgEwDIK+IEf7XMm7x+MM4QAYBkFfkEtnplMfn52Z1r03re/5vH5vEACQhqAvyB3XX67pFVPnPTa9YupcKeVsjzeCXm8QANBLpqC3vcn2YduLtudTjtv2fcnxZ21f1XV8yvZ+21/Nq+FVt3XDrO6+8UrNzkzLal3J333jlefq5fu9EQDAcgxcGWt7StL9kt4v6YikvbZ3RcR3O07bLGld8ucaSQ8kf7d9XNIhSb+SU7srrXsV7L03rb9gQVT76/Z575xeIVu6/dEDumf3YRZRAcgsyxX91ZIWI+KFiHhd0iOStnSds0XSw9HytKQZ26skyfZqSR+Q9IUc211ZnWWVof4LorZumNVT89fp3pvW67Uzb+iVk6cHPgcAumUJ+llJL3V8fSR5LOs5n5H0N5Le6PdNbG+zvWB74fjx4xmaVU1pZZWDqmmGeQ4AtGUJ+rRqv+4NclLPsf2Hkl6OiH2DvklEPBgRcxExt3Jl6m0Pa6FX1Uy/apphngMAbVmC/oikyzq+Xi3paMZzNkr6Y9svqjXkc53tfx26tTXQq2qmXzXNMM8BgLYsQb9X0jrba21fLOlmSbu6ztkl6cNJ9c21kl6NiGMRcWdErI6INcnz9kTEn+XZgaoZppqGChwAoxhYdRMRZ2xvl7Rb0pSkhyLioO1bk+Ofk/SkpBskLUo6KemW8TW52rqrabLsPT/McwCgjf3oAaAG+u1Hzx2mKoY7UQFYLoK+QrrvRNWup5dE2APoiaCfkDyuxPvV0xP0AHoh6Ccgrytx6ukBDIPdKycgr5Wt1NMDGAZBPwF5XYlTTw9gGAT9BOR1JT5oa2MASMMY/QTccf3l543RS8NfiW/dMEuwA1gWgn4CWNkKoEgE/YSM40qcxVMAsiDoK4rFUwCyYjK2orgZCYCsCPqKYvEUgKwI+opi8RSArAj6imLxFICsmIytKEo2AWRF0I/ZOEsgWTwFIAuCfowogQRQBgT9GE1q/3gWTgHoh6Afo0mUQPKpAcAgVN2M0SRKIFk4BWAQgn6MJlECycIpAIMQ9GM0if3jWTgFYBDG6Mds3CWQee51D6CeCPqKY+EUgEEI+hpg4RSAfhijB4CaI+gBoOYYuqkZVskC6EbQ1wirZAGkYeimRlglCyANQV8jrJIFkIagrxFWyQJIQ9CPwc79S9q4Y4/Wzj+hjTv2aOf+pYl8X24vCCANk7E5K3JClFWyANIQ9Dmb1M1GemGVLIBuDN3kjAlRAGVD0OeMCVEAZZMp6G1vsn3Y9qLt+ZTjtn1fcvxZ21clj7/V9rdtP2P7oO2/z7sDZVOWCdGiJoQBlM/AMXrbU5Lul/R+SUck7bW9KyK+23HaZknrkj/XSHog+fs1SddFxC9sr5D0P7b/MyKezrkfpVGGCVFWyALolGUy9mpJixHxgiTZfkTSFkmdQb9F0sMREZKetj1je1VEHJP0i+ScFcmfyK31JVX0hGjRE8IAyiXL0M2spJc6vj6SPJbpHNtTtg9IelnS1yPiW2nfxPY22wu2F44fP56x+UjDhDCATlmC3imPdV+V9zwnIs5GxHpJqyVdbft30r5JRDwYEXMRMbdy5coMzUIvTAgD6JQl6I9Iuqzj69WSji73nIg4Iem/JW1abiOxPGWZEAZQDlmCfq+kdbbX2r5Y0s2SdnWds0vSh5Pqm2slvRoRx2yvtD0jSbanJf2+pP/Lr/lIs3XDrO6+8UrNzkzLkmZnpnX3jVcyPg801MDJ2Ig4Y3u7pN2SpiQ9FBEHbd+aHP+cpCcl3SBpUdJJSbckT18l6UtJ5c5Fkh6LiK/m3w10K3pCGEB5uFUoUy5zc3OxsLBQdDMAoDJs74uIubRjrIwFgJpjU7MG4D6yQLMR9DXHKlkADN3UHPeRBUDQ1xyrZAEQ9DXHKlkABH2Oyrg1MKtkATAZm5OyTnqWYdtkAMUi6HNS5q2BWSULNBtDNzlh0hNAWRH0OWHSE0BZEfQ5qcKkZxkniwGMH2P0OSn7pGdZJ4sBjB9Bn6MyT3qWebIYwHgxdNMQTBYDzUXQNwSTxUBzEfQNUYXJYgDjwRh9Q5R9shjA+BD0DVLmyWIA48PQDQDUHFf0DcXtBYHmIOgbiMVTQLMwdNNA3F4QaBaCvoFYPAU0C0HfQCyeApqFoG8gFk8BzcJkbAOxeApoFoJ+RFUtU2TxFNAcBP0IKFMEUAUE/QjqsMd7VT+RAMiOoB9B1csU+UQCNANVNyOoepkiC6eAZiDoR1D1MsWqfyIBkA1BP4KtG2Z1941XanZmWpY0OzOtu2+8sjLDHlX/RAIgG8boR1TlMsU7rr/8vDF6qVqfSABkQ9A3GAungGYg6Buuyp9IAGRD0OMcauqBeso0GWt7k+3Dthdtz6cct+37kuPP2r4qefwy29+wfcj2Qdsfz7sDyEe7pn7pxCmF3qyp37l/qeimARjRwKC3PSXpfkmbJV0h6YO2r+g6bbOkdcmfbZIeSB4/I+kTEfFbkq6V9LGU56IEqKkH6ivLFf3VkhYj4oWIeF3SI5K2dJ2zRdLD0fK0pBnbqyLiWER8R5Ii4ueSDkliLKCEqKkH6itL0M9Keqnj6yO6MKwHnmN7jaQNkr617FZi7KipB+orS9A75bFYzjm23y7pK5Jui4ifpX4Te5vtBdsLx48fz9As5Knqq3wB9JYl6I9Iuqzj69WSjmY9x/YKtUL+3yLi8V7fJCIejIi5iJhbuXJllrYjR1Vf5QugtyzllXslrbO9VtKSpJslfajrnF2Sttt+RNI1kl6NiGO2LemLkg5FxKdzbDfGgJp6oJ4GBn1EnLG9XdJuSVOSHoqIg7ZvTY5/TtKTkm6QtCjppKRbkqdvlPTnkp6zfSB57G8j4slce4FcUU8P1Isjuofbizc3NxcLCwtFN6ORuveol1pj9QzjAOVme19EzKUdY/dKnId6eqB+CHqch3p6oH7Y62ZIdR3HvnRmWkspoU49PVBdXNEPoc77wlBPD9QPQT+EOo9jU08P1A9DN0Oo+zh2Zz19e4jq9kcP1GqICmgSruiH0JR9Yeo8RAU0CUE/hKaMY9d5iApoEoZuhtCUe63WfYgKaAqCfkhN2BeGUkugHhi6QU9NGaIC6o6gR0/dpZYz0yv01hUX6fZHD2jjjj1MygIVQdCjr60bZvXU/HW696b1eu3MG3rl5GkqcICKIeiRCRU4QHUR9MiEChygugh6ZNKURWJAHRH0yIQKHKC6qKNHJp2LxJZOnNKUfd4Yfd3XFABVxhU9Mtu6Yfbclf3Z5BaUVN8A5UfQY1movgGqh6DHslB9A1QPQY9lofoGqB6CHstC9Q1QPQT9Muzcv6SNO/Zo7fwTjd3rhf1vgOoh6DPibktvYv8boFoI+oyoNrkQPxOgGgj6jKg2uRA/E6AaCPqMqDa5ED8ToBoI+oyoNrlQ2s9Ekk6+foZxeqBE2Osmo6bcEHw52n2/a9dBnTh1+tzjr5w8rTsff+68cwAUx5HsWVImc3NzsbCwUHQzkNHGHXtSbyI+OzOtp+avK6BFQPPY3hcRc2nHGLrByJiUBcqNoMfIek2+hsQiKqAECHqMrNekrMQiKqAMCHqMrHNbhDQsogKKRdAjF+1tEdzjOOP1QHEIeuSq13j9RTbDN0BBCHrkqtd4/dkIxuqBghD0yFV7vH7KFw7iMFYPFCNT0NveZPuw7UXb8ynHbfu+5Piztq/qOPaQ7ZdtP59nw1FeWzfM6o0eC/EYqwcmb2DQ256SdL+kzZKukPRB21d0nbZZ0rrkzzZJD3Qc+2dJm/JoLKqD2nqgPLJc0V8taTEiXoiI1yU9ImlL1zlbJD0cLU9LmrG9SpIi4puSfppnoyeNO0stH7X1QHlkCfpZSS91fH0keWy55/Rle5vtBdsLx48fX85Tx4o7Sw2H2nqgPLIEfVppdPcAbJZz+oqIByNiLiLmVq5cuZynjhV3URoetfVAOWQJ+iOSLuv4erWko0OcU0ls2DU6auuBYmUJ+r2S1tlea/tiSTdL2tV1zi5JH06qb66V9GpEHMu5rYXgLkqjo7YeKNbAoI+IM5K2S9ot6ZCkxyLioO1bbd+anPakpBckLUr6vKS/bD/f9pcl/a+ky20fsf3RnPswVtxZanTU1gPF4sYjGezcv8SdpXKwdv6JnhM3s/xcgZH0u/EItxLMYOuGWQIoB5fOTKfeiUp6s5pJ4vaDQN7YAgET06+2XmIYBxgXgh4TM6i2Xmpd2TM5C+SLoMdEtWvr+4U9lThAvgh6FKLfMA5DOEC+CHoUoj2M08vSiVPsKwTkhKBHYbZumB04Xs8wDjA6gh6FylKJ84nHniHsgRFQR98Di6Qmo/0zvWf34Z419u2tEjrPB5AdK2NTtLcm7ty1cnrFlO6+8UqCZow27tjTM+zbWEELpOu3MpahmxRsTVyMQcM4EuP2wDAI+hRsTVyMfpufdTp1+qxue/QAVTlARgR9CrYmLs7WDbP61J/+7sAre4mreyArgj4FWxMXK8tWCW1U5QCDEfQpOoPGak0AMhE7We2tEj5z0/qBV/fcwAToj6oblF671HVQRY5EVQ6ai6obVNpyru4ZtwcuRNCjMpZTlcO4PfAmhm5QOWkL2tJYUojhHDQDtxJErWTZNkHSufvTcptCNB1B34U9bqqhfR/frFf37eGc9nOBJmHopgN73FTTzv1L+sRjz+hsht9lhnNQV1TdZMQeN9W0nNW0ncM5tz16QBv+4WtM2qL2CPoO7HFTXd2rafvX5bzplZOndfujB7Rm/gn2zkFtEfQd2OOm2tr19i/u+IDuvWn9wDLMNq7yUXcEfQf2uKmP5QzndOMqH3VD0Hdgj5t6GXY4R+IqH/VC1Q0aY+f+Jd2166BOnDo91PPbFTsz0ytkSydOnqYEF6XRr+qGoBe1803TuUlaO7xH9a63rdDf/dFv83uDwhD0fVA732yjXuV3ar9pTNk6G0GtPiaKoO+j1w2pZ2em9dT8dRNpA4o3jqt8ifDH5BD0faydfyL1f2pL+uGOD0ykDSiXPK/y03SHP2P+yANB3wdX9OhlXFf5WbxtxUV6y4opwh+ZEfQp+v1PzBg9unVO2L9zeoVeP3NWJ0+/MdE2tMP/lZOn+TSACxD0XdImYNnsCstV5BV/L72GhXhzqD+CvgvDNchbZ+i3A7Us4Z9Fv08LvElUA0HfYef+Jd326IHUY0zAIk9VD/+sOucT3tnnzWGUx3hjGazxQZ/1IzZX9Bi3tPCfKWjMv6oGzVWM681m1DeqrO0advh45KC3vUnSZyVNSfpCROzoOu7k+A2STkr6i4j4Tpbnphkm6Ef9H4gJWBStDBO+KIdh8mikoLc9Jel7kt4v6YikvZI+GBHf7TjnBkl/pVbQXyPpsxFxTZbnpllu0Ge9nVw/n7lpPSGP0ukO/35XhHUcFmqy5Y4wjHpz8KslLUbEC8l/7BFJWyR1hvUWSQ9H613jadsztldJWpPhuSNLuzPUcszOTBPyKKX2vXGz6PWpNu3NgU8L5ZfnDY+yBP2spJc6vj6i1lX7oHNmMz5XkmR7m6RtkvSe97wnQ7PeNMoPhP3mURfLeVOQlvdpgTeJycvzhkdZgj5tG+/uT4i9zsny3NaDEQ9KelBqDd1kaNc5l85Mp5ZL9kLNPLD8N4ZOaW8S45oIbeIbS94XoFmC/oikyzq+Xi3paMZzLs7w3JHdcf3lmcfo2U4WGN0obxLDyPrpo8lVN/1kCfq9ktbZXitpSdLNkj7Udc4uSduTMfhrJL0aEcdsH8/w3JG1fyC9xiepwQWqbdJvLHUzMOgj4ozt7ZJ2q1Ui+VBEHLR9a3L8c5KeVKviZlGt8spb+j13HB3hFwEA0jViwRQA1F2/8kpuDg4ANUfQA0DNEfQAUHMEPQDUXCknY5OyzB9JukTSTwpuTpGa3H/63lxN7v8off+NiFiZdqCUQd9me6HXLHITNLn/9L2ZfZea3f9x9Z2hGwCoOYIeAGqu7EH/YNENKFiT+0/fm6vJ/R9L30s9Rg8AGF3Zr+gBACMi6AGg5kob9LY32T5se9H2fNHtGTfbL9p+zvYB2wvJY++2/XXb30/+flfR7cyD7Ydsv2z7+Y7HevbV9p3J78Fh29cX0+r89Oj/XbaXktf/QHIf5vax2vTf9mW2v2H7kO2Dtj+ePF77179P38f/2kdE6f6otaXxDyS9V62blzwj6Yqi2zXmPr8o6ZKux/5J0nzy73lJ/1h0O3Pq6/skXSXp+UF9lXRF8vq/RdLa5Pdiqug+jKH/d0n665Rza9V/SaskXZX8+x2Svpf0sfavf5++j/21L+sV/bkbkkfE65LaNxVvmi2SvpT8+0uSthbXlPxExDcl/bTr4V593SLpkYh4LSJ+qNY9D66eRDvHpUf/e6lV/yPiWER8J/n3zyUdUuve0rV//fv0vZfc+l7WoO91s/E6C0lfs70vuVG6JP16RByTWr8kkn6tsNaNX6++Nul3YbvtZ5OhnfbQRW37b3uNpA2SvqWGvf5dfZfG/NqXNegz31S8RjZGxFWSNkv6mO33Fd2gkmjK78IDkn5T0npJxyR9Knm8lv23/XZJX5F0W0T8rN+pKY9Vuv8pfR/7a1/WoM9yQ/JaiYijyd8vS/p3tT6i/dj2KklK/n65uBaOXa++NuJ3ISJ+HBFnI+INSZ/Xmx/Ra9d/2yvUCrp/i4jHk4cb8fqn9X0Sr31Zg/7cDcltX6zWTcV3FdymsbH9y7bf0f63pD+Q9Lxaff5IctpHJP1HMS2ciF593SXpZttvSW4yv07Stwto31i1Qy7xJ2q9/lLN+m/bkr4o6VBEfLrjUO1f/159n8hrX/RMdJ8Z6hvUmpX+gaRPFt2eMff1vWrNrj8j6WC7v5J+VdJ/Sfp+8ve7i25rTv39slofUU+rddXy0X59lfTJ5PfgsKTNRbd/TP3/F0nPSXo2+R98VR37L+n31Bp+eFbSgeTPDU14/fv0feyvPVsgAEDNlXXoBgCQE4IeAGqOoAeAmiPoAaDmCHoAqDmCHgBqjqAHgJr7f66bW2fPBmArAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the posterior probability density function for the paleointensity\n",
    "\n",
    "plt.scatter(test1['Bamb'],test1['normalized'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "417af22f",
   "metadata": {},
   "source": [
    "## Example 2: Model selection, spinning vs non-spinning chondrule in constant ambient field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "c684e697-317f-4469-9c7d-9d75e0c09e29",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bayes Factor in favor of spin\n",
      "3.938989835657311\n",
      "Probability of spin\n",
      "0.7975294476655035\n"
     ]
    }
   ],
   "source": [
    "# Comparison of two models with constant ambient field \n",
    "\n",
    "chond_pints = np.array([37.0, 40.0, 24.5, 11.5, 21.0, 0.0, 0.0])\n",
    "\n",
    "chond_sigmas = np.array([11.0, 13.5, 8.7, 4.7, 5.5, 8.7, 8.7])\n",
    "\n",
    "# run for constant field, no chondrule spin\n",
    "nospin ,Brec,PDF= likelihood_table(chond_pints, chond_sigmas, Btype=0, \n",
    "                                  Bamb_min=1., Bamb_max=100., Bamb_int=1)\n",
    "\n",
    "# version with artificially small \"standard error\" on constant Bamb to simulate Dirac delta\n",
    "#nospin ,Brec,PDF= likelihood_table(chond_pints, chond_sigmas, Btype=0, \n",
    "#                                  Bamb_min=1., Bamb_max=100., Bamb_int=1, PSE=0.1) \n",
    "\n",
    "# run for constant field, yes chondrule spin\n",
    "yesspin ,Brec,PDF= likelihood_table(chond_pints, chond_sigmas, Btype=1, \n",
    "                                  Bamb_min=1., Bamb_max=100., Bamb_int=1)\n",
    "\n",
    "print('Bayes Factor in favor of spin')\n",
    "bayesfactor=yesspin['total']/nospin['total']\n",
    "print(bayesfactor)\n",
    "print('Probability of spin')\n",
    "print(1-(1+bayesfactor)**(-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "4e70faba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fa078a843a0>"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbZklEQVR4nO3df2zc933f8ec7NJWc7UVnx8pineRKBQSmmpWULue6dVB0UTdaShtpQpHKRZYu2KIJjfMLHQMJHTKv2CBjDNrEmGtXdtzVbRpDcDlNS40ygx1gaDF7osxMsqKw0ZQ1IqnOTG3KaXyoKOm9P77fk46n+/E53h3ve597PQDCx++Pu8/HIl/fDz+fz/fzNXdHRETi9bZuF0BERDpLQS8iEjkFvYhI5BT0IiKRU9CLiETupm4XoJo77rjDN23a1O1iiIj0jBMnTvzA3ddV25fJoN+0aRNTU1PdLoaISM8ws7+qtU9dNyIikVPQi4hETkEvIhI5Bb2ISOQU9CIikcvkrJt+cnR6jvHJGeYXi6zP5xgbHWL3cKHbxRKRiFgWV68cGRnxfpheeXR6joMTpyguXbm2zQAHCgp9EWmCmZ1w95Fq+9R100XjkzPLQh6SkAeYWyxycOIUR6fnVr9gIhIVBX0XzS8W6+4vLl1hfHJmlUojIrFS0HfR+nyu4TGNLgYiIo0o6LtobHSI3OBA3WNCLgYiIvVo1k0XlQZaxydnmFssXhuILckNDjA2OtSVsolIPBT0XbZ7uHAt8DXVUkQ6QUGfIeWhLyLSLuqjFxGJnIJeRCRyCnoRkcgp6EVEIqegFxGJnGbddIGmUYrIalLQr7LKFStLi5cBy8JeFwMRaRd13ayyaitWVi5eVroYzC0WcbSSpYi0JijozewBM5sxs7NmdqDK/vea2f80s78zs3/dzLn9ptYiZeXbQy4GIiKhGga9mQ0AjwE7gK3Ag2a2teKw14FPA19cwbl9pdYiZeXbQy4GIiKhQlr09wJn3f2cu18CngV2lR/g7q+5+3Fgqdlz+021FSsrFy8LuRiIiIQKCfoCcL7s+9l0W4jgc81sn5lNmdnUwsJC4Nv3nt3DBQ7t2UYhn8NIHhl4aM+2ZQOtIRcDEZFQIbNurMq20AfNBp/r7oeBw5A8Mzbw/XtSo8XLypcv1qwbEWlVSNDPAhvLvt8AzAe+fyvn9jWtZCki7RLSdXMc2GJmm81sDbAXOBb4/q2cKyIibdCwRe/ul83sIWASGACedvfTZrY/3f+Emb0HmALeCVw1s88CW939zWrndqguIiJShblnrzt8ZGTEp6amul0MEZGeYWYn3H2k2j7dGSsiEjkFvYhI5BT0IiKRU9CLiEROQS8iEjkFvYhI5BT0IiKRU9CLiEROQS8iEjkFvYhI5BT0IiKRU9CLiEROQS8iEjkFvYhI5EKeMCVtcHR6To8GFJGuUNCvgqPTcxycOEVx6QoAc4tFDk6cAggKe10kRKQV6rpZBeOTM9dCvqS4dIXxyZmG55YuEnOLRZzrF4mj03MdKq2IxEZBvwrmF4tNbS/XykVCRAQU9KtifT7X1PZyrVwkRERAQb8qxkaHyA0OLNuWGxxgbHSo4bmtXCREREBBvyp2Dxc4tGcbhXwOAwr5HIf2bAsaUG3lIiEiApp1s2p2DxdWNFOmdI5m3YjISinoe8BKLxIiIqCuGxGR6CnoRUQip6AXEYmcgl5EJHJBQW9mD5jZjJmdNbMDVfabmT2a7j9pZveU7fucmZ02s1fN7Gtm9o52VkBEROprGPRmNgA8BuwAtgIPmtnWisN2AFvSr33A4+m5BeDTwIi73w0MAHvbVnoREWkopEV/L3DW3c+5+yXgWWBXxTG7gGc88RKQN7M70303ATkzuwm4GZhvU9lFRCRASNAXgPNl38+m2xoe4+5zwBeB7wMXgIvu/o1qH2Jm+8xsysymFhYWQssvIiINhAS9VdnmIceY2W0krf3NwHrgFjP7aLUPcffD7j7i7iPr1q0LKJaIiIQICfpZYGPZ9xu4sful1jG/AHzP3RfcfQmYAH525cUVEZFmhQT9cWCLmW02szUkg6nHKo45BnwsnX1zH0kXzQWSLpv7zOxmMzNgO3CmjeUXEZEGGq514+6XzewhYJJk1szT7n7azPan+58Angd2AmeBt4CPp/teNrPngFeAy8A0cLgTFRERkerMvbK7vftGRkZ8amqq28UQEekZZnbC3Ueq7dOdsSIikVPQi4hETkEvIhI5Bb2ISOQU9CIikVPQi4hETkEvIhI5Bb2ISOQU9CIikVPQi4hEruFaN7JyR6fnGJ+cYX6xyPp8jrHRIXYPVy7lLyLSWQr6Djk6PcfBiVMUl64AMLdY5ODEKYCWwl4XDxFplrpuOmR8cuZayJcUl64wPjmz4vcsXTzmFos41y8eR6fnWiytiMRMQd8h84vFpraH6MTFQ0Tip6DvkPX5XFPbQ3Ti4iEi8VPQd8jY6BC5wYFl23KDA4yNDq34PTtx8RCR+CnoO2T3cIFDe7ZRyOcwoJDPcWjPtpYGTjtx8RCR+GnWTQftHi60dUZM6b0060ZEmqGg7zHtvniISPzUdSMiEjkFvYhI5BT0IiKRU9CLiEROQS8iEjkFvYhI5BT0IiKRU9CLiEQuKOjN7AEzmzGzs2Z2oMp+M7NH0/0nzeyesn15M3vOzL5jZmfM7GfaWQEREamvYdCb2QDwGLAD2Ao8aGZbKw7bAWxJv/YBj5ft+zLwZ+7+XuD9wJk2lFtERAKFtOjvBc66+zl3vwQ8C+yqOGYX8IwnXgLyZnanmb0T+DngKwDufsndF9tXfBERaSQk6AvA+bLvZ9NtIcf8OLAA/L6ZTZvZU2Z2S7UPMbN9ZjZlZlMLCwvBFRARkfpCgt6qbPPAY24C7gEed/dh4EfADX38AO5+2N1H3H1k3bp1AcUSEZEQIUE/C2ws+34DMB94zCww6+4vp9ufIwl+ERFZJSFBfxzYYmabzWwNsBc4VnHMMeBj6eyb+4CL7n7B3f8aOG9mpSdjbAe+3a7Ci4hIYw3Xo3f3y2b2EDAJDABPu/tpM9uf7n8CeB7YCZwF3gI+XvYWnwK+ml4kzlXsExGRDjP3yu727hsZGfGpqaluF0NEpGeY2Ql3H6m2T3fGiohETkEvIhI5Bb2ISOQU9CIikVPQi4hETkEvIhI5Bb2ISOQa3jAl2XV0eo7xyRnmF4usz+cYGx1i93DlenMi0u8U9G22WuF7dHqOgxOnKC5dAWBuscjBiVMACnsRWUZdN21UCt+5xSLO9fA9Oj3X9s8an5y5FvIlxaUrjE/OtP2zRKS3KejbaDXDd36x2NR2EelfCvo2Ws3wXZ/PNbVdRPqXgr6NVjN8x0aHyA0OLNuWGxxgbHSoxhki0q8U9G20muG7e7jAoT3bKORzGFDI5zi0Z5sGYkXkBpp100alkF2tKY+7hwsKdhFpSEHfZgpfEckadd2IiEROQS8iEjkFvYhI5BT0IiKRU9CLiEROQS8iEjkFvYhI5BT0IiKRU9CLiEROQS8iEjkFvYhI5IKC3sweMLMZMztrZgeq7DczezTdf9LM7qnYP2Bm02b29XYVXEREwjQMejMbAB4DdgBbgQfNbGvFYTuALenXPuDxiv2fAc60XFoREWlaSIv+XuCsu59z90vAs8CuimN2Ac944iUgb2Z3ApjZBuBDwFNtLLeIiAQKCfoCcL7s+9l0W+gxXwI+D1xdWRFFRKQVIevRW5VtHnKMmf0i8Jq7nzCzn6/7IWb7SLp9uOuuuwKKJZl08gi88FtwcRZytyXbim8sf712A2z/ArzvI90rp0gfCQn6WWBj2fcbgPnAY34Z+LCZ7QTeAbzTzP7I3T9a+SHufhg4DDAyMlJ5IZEsuxbu50mu+ek/X/H168eUv754Hib2wcQnYO1Ghb5Ih4V03RwHtpjZZjNbA+wFjlUccwz4WDr75j7gortfcPeD7r7B3Tel571YLeSlh508Av/t02nIw41/7NWSHnfxfHL+ySOdKJ2IEBD07n4ZeAiYJJk5c8TdT5vZfjPbnx72PHAOOAs8Cfx6h8orWXHyCPzO3UmrfKnY2nstFZP3+Z27FfgiHWDu2eslGRkZ8ampqW4XQ2opteJbDfhqBnPwS4+qK0ekSWZ2wt1Hqu3Tw8EjcXR6jvHJGeYXi6zP5xgbHercQ8pf+K3OhDwk7/tf9id9+Bq0FWkLBX0brGrI1vj8gxOnKC5dAWBuscjBiVMA7S3HskHXetIB2dztybfLZt28zrIB22o8qce1/ntQ2Iu0QF03LaoMWYDc4ACH9mxbtbC//5EXmVu8sYVdyOf4iwMfbM+HhHbXhMyiCb5gNPm+In2sXteNFjVr0fjkzLKQByguXWF8cmbVyjBfJeTrbV+RRt01gznY8yR87tXGYfy+jyTH7XkyOS+EZueIrJiCvkWrErINrM9XD8ta25tSml1Tr/W9duPKBlDf95HkvLUbAQMbqH/8UjG54IhIUxT0LepoyAYaGx0iN7g8JHODA4yNDrX2xjfMka9i7cawVnwtpdb9w4vwT59o3MK/eF7TMEWapKBvUcdCtgm7hwsc2rONQj6HkfTNt2WMIKS7ZvsXWvuMcsta+HWoG0ekKRqMbYNuz7rpmIfz1Jwd0+nB0ZDB39JfEyKiefSdtnu4EEewV1q7oXq3zWoEbOkCUm92TqkbR7NxROpS143caNkAbMXCpO3urqmn1H9frytH3TgiDSnoZbmqi5SlYb/S2TWt2v6F+oO0mo0jUpeCXparOgDrrc+uaUXIIK1m44jUpKCX5S7ONrd9tagbR2TFFPSy3NoNzW1fberGEWmagl4SWRmAbUTdOCJNU9BLNgdg61E3jkhTFPSSzQHYEOrGEQmioJfsDsA2om4ckSAKesn+AGw96sYRaUhB3896ZQA2hLpxRGrSWjf96oZFw0oDsN6bT3MKWhsn411RIh2iFn2/6tUB2HoaduO4+uulLyno+1WvDsCGqNeNo/566UMK+ggdnZ7j/kdeZPOBP+X+R17k6PTcjQf18gBsI41m46i/XvqMgj4yR6fnODhxirnFIg7MLRY5OHHqxrCv1urttQHYekrdOJWDzCWadil9RIOxK5DlJ0qNT85QXLqybFtx6QrjkzNJGU8eSQcsZyF3G9yUg+IbSUu+1wZgQ9R6eApc78aB+OotUkYt+iYFt5i7ZH6x+qP35heLFUsdOBRfh8tF2HO4dwdgG9G0SxEFfbPqtZizYH2+eqitz+eqz7SJPeh096xIWNCb2QNmNmNmZ83sQJX9ZmaPpvtPmtk96faNZvZNMztjZqfN7DPtrsBqq9tizoCx0SFygwPLtuUGBxgbHYp7pk09untW+lzDoDezAeAxYAewFXjQzLZWHLYD2JJ+7QMeT7dfBn7D3X8CuA/4ZJVze0rdFnMG7B4ucGjPNgr5HAYU8jkO7dmW9M/HPNMmhLpxpE+FDMbeC5x193MAZvYssAv4dtkxu4Bn3N2Bl8wsb2Z3uvsF4AKAu//QzM4AhYpze8rY6BAHJ04t67651mLOiN3DheqDw9u/UHE3LHHNtGkk6O7ZtBsnxoFp6VshXTcFoPy3Yjbd1tQxZrYJGAZervYhZrbPzKbMbGphYSGgWN1Rt8WcVaU1bSb2JbNscrcDls215jtN3TjSh0Ja9NUmInszx5jZrcCfAJ919zerfYi7HwYOA4yMjFS+f6bUbDFnUeWaNsXXk1b8nsP9FfCVqv11U67UjdPP/48kGiEt+lmgvPmzAZgPPcbMBklC/qvuPrHyosqK9ONMmxCajSN9JCTojwNbzGyzma0B9gLHKo45BnwsnX1zH3DR3S+YmQFfAc64+2+3teQSpl9n2oRQN470iYZB7+6XgYeASeAMcMTdT5vZfjPbnx72PHAOOAs8Cfx6uv1+4J8BHzSzb6VfO9tdCamj32fahAiZjTPxCbXupWcFLYHg7s+ThHn5tifKXjvwySrn/Tk1FxuRVdHvM21ChMzGAS2ZID1Ld8bGqDTL5uF8El7v/9W0e6JPZ9qECOnGAY1vSE/SomaxqZxlc/E8/O8/VriHajQbB5L/pw/n410ITqKjFn1sNMumNSGzcQBwDdRKz1DQx6bGbJqrF2frP4hErit14+x5sv4gLWigVnqCgj42NWbTzF99VyaXVc60Za37BnMK1LqXDFPQBwp6PF83lQZgL56nMpTe8jX8x8vX+5GztKxy5pVa9w8vhg3UqnUvGaSgD5D1h40sf6AIJKtPJGE/e/UODiz9S45d/cCyU7KyrHJPaTTfvkSte8kYBX2ArD9spOoALA5rN/IrNz95Q8hDdpZV7inBA7WodS+ZoqAPkPWHjdRb5qDug0ikec0M1IJa95IJCvoAWX/YSL1lDnpyWeVeoNa99BDdMBUgsw8bOXmk7LZ9Y9nq0WXLHPTUssq95H0fSb4qb1KrRUsoSJeoRR8gk63iOgOwWuZglal1LxlnyXpk2TIyMuJTU1PdLka2XZtKWWHtxqQPWbojtHUPXPsrbO1GLaUgLTOzE+4+Um2fWvS9SuvMZ1MzrftSV5sGbKXDFPS9pnRj1A1Pc0w1WGc+8zd+xaDZmTmgLh3pKA3G1nF0eo7xyRnmF4usz+cYGx3KRr98rW6BBuvMl278Kg0ql278AjRY2wmh69yXu3g+eYj7xCfUpSNtoxZ9DZm8G7bqjVGpgAHYzN/4FaOVtO7Lu3Qm9sHDa9XSl5Yo6GvIZCjW7H+3JEwatPwyf+NXzG7ouw998JpCX1qnoK8hU6HYYr98SeZv/IrdtQXSLsKew4EDtuUU+rIyCvoaMhOKN8yXr9DE81+1HEKGrKhLp5xCX8JpHn0NlQOXkITiqt8oVWu+PKxosK58gHltbhAzWHxrKRuDzf2q3h3OTSubm7/ln8B3v5F0+emxh9GrN49eQV8hM0G47Je/GkvWSF+hzFzIZLm2hn4l3aAVMwV9oMyEX8jdlS3eAXv/Iy8yV2W8oZDP8RcHPrji95U2Wo3Qz92efFt8Q63+Hlcv6DWPvky9mTarEvQNW/GpJvrla8nUYLNUV1o0DToQ+un5xdevbyqfw19+Acjddv21LgY9SS16rnfXVGvhQvJr9b1HPtSZD2/2F7hNf3bXatFD0qpXf32GdbSlH6LKXwPlFwNdGLpCXTd1VOuuqdT27oyV/qK2ccGyRvUulUqhn3FdD/0QNbqJygeLa10odNEIpqCvolErvqSlPvprv4TlP8ivs6JfyMFc25ceDv1/oNDvEeU/b8uCNMsXgWY18ddEu16v5KLUhQtUy0FvZg8AXwYGgKfc/ZGK/Zbu3wm8Bfxzd38l5NxqVhL0x4/9HhtfGefdvsBFuxUw1voPa75+w2/FDPL8bc3Xb9rfIzc4wNuXLq7gh2CFgV5Nh2dJbD7wp8GlLNUoXzYjaW2bXq/P5/hH713HN7+zcMOspyy8znr56pV17zte4lP8Me/xH4DpBppmlT3tYUWuAuawGJBNr9k6zt8zxj/88L9q6jNaCnozGwD+EvjHwCxwHHjQ3b9ddsxO4FMkQf/TwJfd/adDzq2m2aA/fuz3uPvEvyFnl4LP6QkdaMVXU6+/XuLz4bf9OZ+/6Qjr7W94w2+51rABeFsraSZtU/Q1vPpT/76psG911s29wFl3P5e+2bPALqA8rHcBz3hy1XjJzPJmdiewKeDclm18ZTyikF/9uc7VHpUo8Tp29QMcu/SBG7ZXvwD8SBeDLsjZJTa+Mg5NtuprCQn6AlA+32+WpNXe6JhC4LkAmNk+YB/AXXfdFVCs697tC639XdV13b2RpdTnXuqvj6U3V5pT6wJQrvHFIPzC4A7W07+3nfVu/0Hb3isk6Kv9U1TmQK1jQs5NNrofBg5D0nUTUK5rXrN1vIeFZk7JgGzdpVj+APHyQVqFvpQLuRiUq3VhmPd38cLVn2T7277VlotGp3XjovSa3cF72vReIUE/C5Qvs7cBmA88Zk3AuS07f88YazPbR997dyAq9KVdGl0Y/m0T79XcXxPte938Ran1C1TR13D+p8baFvQhg7E3kQyobgfmSAZUf9XdT5cd8yHgIa4Pxj7q7veGnFtNa7NufsBFu4VkBPtva75+025tPKOm1alXGQ/0ZtVaB6hfZrVkvXy9VNasl6+dZb0+4+lvgrLpNbtj9WfdpG+wE/gSyRTJp939P5jZfgB3fyKdXvmfgAdIpld+3N2nap3b6POysHqliEgv0Q1TIiKRqxf0um9CRCRyCnoRkcgp6EVEIqegFxGJXCYHY81sAfirJk65A2jfbWS9oR/rDP1Z736sM/RnvVup84+5+7pqOzIZ9M0ys6lao82x6sc6Q3/Wux/rDP1Z707VWV03IiKRU9CLiEQulqA/3O0CdEE/1hn6s979WGfoz3p3pM5R9NGLiEhtsbToRUSkBgW9iEjkejrozewBM5sxs7NmdqDb5ekUM9toZt80szNmdtrMPpNuv93M/ruZfTf9723dLmu7mdmAmU2b2dfT7/uhznkze87MvpP+m/9M7PU2s8+lP9uvmtnXzOwdMdbZzJ42s9fM7NWybTXraWYH03ybMbPRlX5uzwZ9+uDxx4AdwFbgQTPb2t1Sdcxl4Dfc/SeA+4BPpnU9ALzg7luAF9LvY/MZ4EzZ9/1Q5y8Df+bu7wXeT1L/aOttZgXg08CIu99NsqT5XuKs838mWc69XNV6pr/je4F/kJ7zu2nuNa1ng56yh5a7+yWg9ODx6Lj7BXd/JX39Q5Jf/AJJff8gPewPgN1dKWCHmNkG4EPAU2WbY6/zO4GfA74C4O6X3H2RyOtN8rS7XPqwoptJnkQXXZ3d/X8Ar1dsrlXPXcCz7v537v494CxJ7jWtl4O+1gPJo2Zmm4Bh4GXg77v7BUguBsC7u1i0TvgS8Hngatm22Ov848AC8Ptpl9VTZnYLEdfb3eeALwLfBy4AF939G0Rc5wq16tm2jOvloA9+8HgszOxW4E+Az7r7m90uTyeZ2S8Cr7n7iW6XZZXdBNwDPO7uw8CPiKPLoqa0T3oXsBlYD9xiZh/tbqkyoW0Z18tBH/LQ8miY2SBJyH/V3SfSzf/PzO5M998JvNat8nXA/cCHzez/knTLfdDM/oi46wzJz/Wsu7+cfv8cSfDHXO9fAL7n7gvuvgRMAD9L3HUuV6uebcu4Xg7648AWM9tsZmtIBi2OdblMHZE+k/crwBl3/+2yXceAX0tf/xrwX1e7bJ3i7gfdfYO7byL5t33R3T9KxHUGcPe/Bs6b2VC6aTvwbeKu9/eB+8zs5vRnfTvJOFTMdS5Xq57HgL1m9nYz2wxsAf7Xij7B3Xv2C9gJ/CXwf4Df7HZ5OljPD5D8yXYS+Fb6tRN4F8ko/XfT/97e7bJ2qP4/D3w9fR19nYGfBKbSf++jwG2x1xv4d8B3gFeBPwTeHmOdga+RjEMskbTY/0W9egK/mebbDLBjpZ+rJRBERCLXy103IiISQEEvIhI5Bb2ISOQU9CIikVPQi4hETkEvIhI5Bb2ISOT+P51HL0PKwrQaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(nospin['Bamb'],nospin['normalized'])\n",
    "plt.scatter(yesspin['Bamb'],yesspin['normalized'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "412ad844",
   "metadata": {},
   "source": [
    "## Example 3: Model selection, spinning vs non-spinning chondrule in sinusoidal ambient field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "ac91e677",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bayes Factor in favor of spin\n",
      "1.8418803359682752\n",
      "Probability of spin\n",
      "0.6481203000198517\n"
     ]
    }
   ],
   "source": [
    "# Comparison of two models with sinusoidal ambient field\n",
    "\n",
    "chond_pints = np.array([37.0, 40.0, 24.5, 11.5, 21.0, 0.0, 0.0])\n",
    "\n",
    "chond_sigmas = np.array([11.0, 13.5, 8.7, 4.7, 5.5, 8.7, 8.7])\n",
    "\n",
    "# run for constant field, no chondrule spin\n",
    "nospin ,Brec,PDF= likelihood_table(chond_pints, chond_sigmas, Btype=2, \n",
    "                                  Bamb_min=10., Bamb_max=250., Bamb_int=1.)\n",
    "# run for constant field, yes chondrule spin\n",
    "yesspin ,Brec,PDF= likelihood_table(chond_pints, chond_sigmas, Btype=3, \n",
    "                                  Bamb_min=10., Bamb_max=250., Bamb_int=1.)\n",
    "\n",
    "print('Bayes Factor in favor of spin')\n",
    "bayesfactor=yesspin['total']/nospin['total']\n",
    "print(bayesfactor)\n",
    "print('Probability of spin')\n",
    "print(1-(1+bayesfactor)**(-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "8dfc2864-4c89-4bce-b293-9edbabb58e4c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fa0b9f49700>"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeLElEQVR4nO3df5Ac5X3n8fdXy8oefpwWghxLIxEpVypcOJK9eA1UiUqV8dlIWKC1UqXDuRyOTSIrgYoxPg5x9hmFSoIwFcWhQsEJQ8Wc7cBeAZs1xhEUcOXCFdlaIbEgsGIFbEsrVZAPS+bHxvr1vT9mRhrNds/0zHTP9HR/XlVbu9v9DPs0s/r2s9/+Ps9j7o6IiGTXjG53QEREkqVALyKScQr0IiIZp0AvIpJxCvQiIhl3Wrc7EOTcc8/1BQsWdLsbIiI9Y9u2bb9w99lB51IZ6BcsWMD4+Hi3uyEi0jPM7Gdh55S6ERHJOAV6EZGMU6AXEck4BXoRkYxToBcRybhUVt2kxej2Se7cvIt9B6eYO1DgpsvPZ3iw2O1uiYg0RSP6EKPbJ7nl0ReZPDiFA5MHp7jh4R0M3vYko9snu909EZHIFOhD3Ll5F1NHjk07/st3jnDLoy8q2ItIz1CgDzF5cCr03NSRY9y5eVcHeyMi0joF+gCj2yexBm321bkRiIikiQJ9gDs376LRvltzBwod6YuISLsU6AM0Gq0b8JH3Ba4dJCKSOpECvZktM7NdZrbbzNYFnDczu6t8fsLMLqw691Mze9HMdphZT6xU1mi07sAj2yb1QFZEekLDQG9mfcDdwHLgAuBTZnZBTbPlwKLyxxrgnprzH3H3D7r7UPtdTt5Nl59Pob+vbhs9kBWRXhFlRH8RsNvdX3X3w8BDwMqaNiuBB71kCzBgZnNi7mvHDA8WuX3VYooDhboPZfVAVkR6QZRAXwT2VH2/t3wsahsHnjSzbWa2JuyHmNkaMxs3s/EDBw5E6FayhgeL/GDdZby24RMUQ1I5eiArIr0gyhIIQYPa2qKUem2Wuvs+M3sP8JSZ/djdvz+tsfsmYBPA0NBQo6KXRNUuffCR983mkW2Tp0ygKvT3cdPl53exlyIi0UQZ0e8F5ld9Pw/YF7WNu1c+vw48RikVlFpBSx88sm2S3/tQ8UQqpzhQ4PZVi7XujYj0hCiBfiuwyMwWmtlM4GpgrKbNGHBNufrmEuCQu+83szPM7CwAMzsD+DjwUoz9j13Q0gdTR47x7I8PcNPl5zN3oMC+g1PcuXmXqm5EpCc0TN24+1Ezux7YDPQBD7j7TjNbWz5/L/AEcAWwG3gH+Ez55b8JPGZmlZ/1bXf/p9ivIkZhD1gnD05xy6MvnrgJVL4HNLIXkVSLtEyxuz9BKZhXH7u36msHrgt43avAB9rsY0fNHSgErnPTZxY40r9z8y4FehFJNc2MrRFUQ1/o7+OYBz8fVomliKSdAn2N2hr6yoNXlViKSK/SDlMBhgeLgemY6hw9qMRSRHqDAn1ElcCvrQVFpNco0DchbKQvIpJmytGLiGScAr2ISMYpdVOjdp2b2jx8o/MiImmjQF+lss5N2OzXRudFRNJIqZsqYevcVDYYaXReRCSNFOirhM1yrRxvdF5EJI0U6KuEzXKtHG90XkQkjRToq4Stc1OZ/drovIhIGulhbJVGs181O1ZEepF5yKqM3TQ0NOTj4+Pd7oaISM8ws23uPhR0TqkbEZGMU6AXEck4BXoRkYzTw9gWaSkEEekVCvQt0FIIItJLlLppgZZCEJFeokDfAi2FICK9RKmbsmZy7nMHCkwGBHUthSAiaaQRPSdz7pMHp3BO5txHt08GttdSCCLSSxToaT7nPjxY5PZViykOFDCgOFDg9lWL9SBWRFJJqRtay7lro3AR6RUa0aPlh0Uk2xToUc5dRLJNqRu0/LCIZFukQG9my4C/BfqAr7v7hprzVj5/BfAO8Ifu/nzV+T5gHJh09xUx9T1WyrmLSFY1DPTlIH038DFgL7DVzMbc/eWqZsuBReWPi4F7yp8rPg+8AvyHmPqdClrvRkR6QZQc/UXAbnd/1d0PAw8BK2varAQe9JItwICZzQEws3nAJ4Cvx9jvrmu29l5EpFuiBPoisKfq+73lY1HbfA3478Dxej/EzNaY2biZjR84cCBCt7pL692ISK+IEugt4Fjt/oOBbcxsBfC6u29r9EPcfZO7D7n70OzZsyN0q7u03o2I9IoogX4vML/q+3nAvohtlgJXmdlPKaV8LjOzb7bc2xRR7b2I9IoogX4rsMjMFprZTOBqYKymzRhwjZVcAhxy9/3ufou7z3P3BeXXPePufxDnBXSLau9FpFc0rLpx96Nmdj2wmVJ55QPuvtPM1pbP3ws8Qam0cjel8srPJNfldFDtvYj0CnOvTbd339DQkI+Pj3e7GyIiPcPMtrn7UNA5LYEgIpJxCvQiIhmntW7QDFcRybbcB/rKDNfK5KfKDFdAwV5EMiH3qRvNcBWRrMv9iL7dGa5K+4hI2uV+RN/ODFctbCYivSD3gb6dGa5K+4hIL8h96qadGa5a2ExEekHuAz20vrvU3IECkwFBXQubiUia5D510w4tbCYivUAj+jZoYTMR6QUK9G3SpuIiknZK3YiIZJwCvYhIxinQi4hknAK9iEjG6WFsDLTejYikmQJ9m7TMsYiknVI3bdJ6NyKSdgr0bdJ6NyKSdgr0bWpnmWMRkU5QoG+T1rsRkbTL/cPYditmtN6NiKRdrgN9XBUzWu9GRNIs16kbVcyISB7kOtCrYkZE8iDXgV4VMyKSB7kO9KqYEZE8iBTozWyZme0ys91mti7gvJnZXeXzE2Z2Yfn4u83sR2b2gpntNLM/j/sC2jE8WOT2VYspDhQwoDhQ4PZVi1t+sDq6fZKlG55h4brvsnTDM4xun4y3wyIiLWhYdWNmfcDdwMeAvcBWMxtz95ermi0HFpU/LgbuKX/+NXCZu79lZv3Ac2b2PXffEvN1tCyuihmteSMiaRVlRH8RsNvdX3X3w8BDwMqaNiuBB71kCzBgZnPK379VbtNf/vC4Op8mquARkbSKEuiLwJ6q7/eWj0VqY2Z9ZrYDeB14yt1/GPRDzGyNmY2b2fiBAwcidj89VMEjImkVJdBbwLHaUXloG3c/5u4fBOYBF5nZ7wT9EHff5O5D7j40e/bsCN1KF1XwiEhaRQn0e4H5Vd/PA/Y128bdDwL/F1jWbCd7gSp4RCStogT6rcAiM1toZjOBq4GxmjZjwDXl6ptLgEPuvt/MZpvZAICZFYD/BPw4vu6nR9wVPCIicWlYdePuR83semAz0Ac84O47zWxt+fy9wBPAFcBu4B3gM+WXzwG+Ua7cmQGMuPvj8V9GOmjNGxFJI3NPXxHM0NCQj4+Pd7sbIiI9w8y2uftQ0Llcz4wVEckDBXoRkYxToBcRybhcbzyShHZ3rBIRiZsCfYy03o2IpJFSNzHSejcikka5HdEnkWLRejcikka5HNFXUiyTB6dwTqZY2l0/XuvdiEga5TLQJ5Vi0Xo3IpJGuUzdJJViqaR+VHUjImmSy0A/d6DAZEBQjyPFovVuRCRtcpm6UYpFRPIklyN6pVhEJE9yGegh2RSLZseKSJrkNtAnRbNjRSRtcpmjT5Jmx4pI2ijQx0yzY0UkbRToY6bZsSKSNgr0MVPppoikjR7GxkylmyKSNgr0CdDsWBFJE6VuREQyToFeRCTjlLpJkGbIikgaKNAnRDNkRSQtlLpJiGbIikhaaESfkKZmyE6MwPduhqk3Th6zGeDHYdZ8+OhXYMnqhHoqIlmXy0Dfidx5pM1NggJ8hR8vfT60Bx79Y3j8BljxNQV8EWla7lI3SW0MXqvhDNnHb4RH1wQH+SCH34Z/vK50cxARaUKkQG9my8xsl5ntNrN1AefNzO4qn58wswvLx+eb2bNm9oqZ7TSzz8d9Ac3qVO58eLDI7asWUxwoYEBxoMDtqxaX/nKYGIHxBwBv7j967HDpLwARkSY0TN2YWR9wN/AxYC+w1czG3P3lqmbLgUXlj4uBe8qfjwJfdPfnzewsYJuZPVXz2o7q5OqSoTNkv3czTQf5iqk34I6FsPwOpXFEJJIoI/qLgN3u/qq7HwYeAlbWtFkJPOglW4ABM5vj7vvd/XkAd38TeAXoam1h11eXnBiJnq4JM/VGKe3z+I3x9ElEMi1KoC8Ce6q+38v0YN2wjZktAAaBHwb9EDNbY2bjZjZ+4MCBCN1qTddXl4wt9eKl9I9y9iLSQJSqGws4Vpt3qNvGzM4EHgFucPdfBf0Qd98EbAIYGhpqMa/RWKdXl6yu8Pn0mT/i1qNvBP7PAmDoWlixsfT1xAg89rmT1TeBvHTjUApHROqIEuj3AvOrvp8H7Ivaxsz6KQX5b7n7o613NT6dWl2yenbsVTOe48tH7sXConzhnJNBHk4G7+/8GRyp8/xg6o3STUHBXkRCREndbAUWmdlCM5sJXA2M1bQZA64pV99cAhxy9/1mZsD9wCvuvpGcqVT4XDXjOTb0f53TrM7ofPkd048tWQ1X3lW6CdSjShwRqaNhoHf3o8D1wGZKD1NH3H2nma01s7XlZk8ArwK7gfuAPy0fXwr8V+AyM9tR/rgi7otIq0olz/r+BzndDoc3LJwTPiJfshpufq2U1glTGdWLiASINDPW3Z+gFMyrj91b9bUD1wW87jmC8/e5MHegwId+9RRn81Z4o/5C8Gi+1oqNsPOx8Iqdx8r3XKVwRKRG7mbGdtJNl5/Pzf0j4Xl56yulZqIG53o3BD+mkksRCaRAn6DhwSJz7RfhDT55b3Mj8CWrG+TrVXIpItMp0CdpYgQLy1zVy8vXs/yOUronlOvhrIicQoE+SU/fRtBSBw7R8vJBKpU41hfeRg9nRaSKAn2SDu0NPOwOo8eWtv7fXbK6lPap95xbo3oRKVOgT1Lh7MDD+/zc9lfLXLIahj4bfl6jehEpU6BPysQI/PrNaYcP+2l89ejqeFbLXLGx/sPZp29r/2eISM9ToE/K07fB8SPTDr/p72bs+KXxrZZZL9d/aI9G9SKSv0A/un2SpRueYeG677J0wzOx7yx1Qkh+/mx7O97VMhuVXH7nzxTsRXIuV4G+U9sIMjFS2tw7wOt27smdpuJSr+TyyJQezIrkXK4CfUe2EZwYKY2i/dj0c/0F3rvqr+JfObNSchlGD2ZFci1Xgb4j2wg+fVvgssLHbQbr/XMs/PYZyaSMlqyGWfPDz2tUL5JbuQr0HdlG8NCe4OPHnb9/66JkU0Yf/Ur4OY3qRXIrV4E+8W0EJ0YIm8S0z3/jlO9jTxlB4wezj61VsBfJoVwF+uHBIrevWkxxoIABxYFCvA9GQ5Y8OO7w1aPT17WJNWVU0WiFS1XhiOROpPXosyTRbQRDSioxGDt+6bTDsaaMKpasLuXjw9atr1ThaN16kdzI1Yg+cbPmBR7+98KcZFNGtRqtcKl8vUiuKNDHadHHmZaj7y9w+vLbkk0Z1YqywqWqcERyQ4E+LhMj8MK3OTVHb/CB34clqxkeLPKDdZfxN//5gwB84eEdyc7MPbHCZQiN6kVyQ4E+LoH18w4/efLEdx2bmVvRqApHo3qRXFCgj0vYg9iq4x2ZmVurXhXO1Btwx0KN7EUyToE+DnXWtql+QNuRmbm1Go3qp95QyaVIxinQt6vB2jbVs1U7MjM3SKNtC7XwmUimKdC3K2RtG6yvVPlSVa+e+MzcMI1G9aCHsyIZpkDfrrDcvB+fNimpMjN3oNB/4ti7+zv0FjSqrQeN6kUyKleBPpFNR0L2hQ2bPAXw66PHT3z9y3eOJFt5U1GprW+Ur9eoXiRzchPoEyltDNkXlr6ZoStJdqXypmLJarj5NZVciuRMbgJ9IgE2ZF9YZp4ZupZMVypvajUqudSoXiRTchPoEwmwYfn5qV+GvqRrlTfVtJyxSK7kJtAnEmDD8vB18vNdq7yp1Wg540fXwOM3dq4/IpKYSIHezJaZ2S4z221m6wLOm5ndVT4/YWYXVp17wMxeN7OX4ux4sxIJsCGLmNXb6amrlTfVGpZcOow/oJG9SAY0jDBm1gfcDSwHLgA+ZWYX1DRbDiwqf6wB7qk69/fAsjg6247YNx1psIhZI12pvKnVsOTS9XBWJAOibDxyEbDb3V8FMLOHgJXAy1VtVgIPursDW8xswMzmuPt+d/++mS2Iu+OtiHXTkQiLmIWp92A4saWLg1RuSI+tDZ7ZC6WHs4/fCCs2dq5fIhKrKDmDIlC94/Xe8rFm29RlZmvMbNzMxg8cONDMS7sjwiJmYVJReVNxYjnj4L1uAaVwRHpclEAfFAFqN0aN0qYud9/k7kPuPjR79uxmXtodLUyUqkhF5U21Jath6LN1GrgqcUR6WJRAvxeYX/X9PGBfC22yo4WJUtWCHgwb8JH3dfEGt2Jj/YezqsQR6VlRAv1WYJGZLTSzmcDVwFhNmzHgmnL1zSXAIXffH3Nf06OFiVLVhgeL/N6Hiqf8GeTAI9smO/9AttryO6ibwsFh/H4Fe5Ee0zDQu/tR4HpgM/AKMOLuO81srZmtLTd7AngV2A3cB/xp5fVm9g/APwPnm9leM7s25mvovBYmStV69scHpuW2OrYUQpgTKZx6wR7l7EV6TJSqG9z9CUrBvPrYvVVfO3BdyGs/1U4H4zC6fZI7N+9i38Ep5g4UuOny89srq7QZwVUqEfLzFal6IFttxUY475L6lTiVnD1E+gtGRLor8zNjY13MrIlNRhoJe/A6w6y76RuIVomjnL1Iz8h8oI91MbMmNhlpJOiBLMAx9+5MnqrVsBIHlLMX6Q2ZD/Sxpkia2GSkkcpM3T6bPmrueq6+YsVGGLqWxjl7BXuRNMt8oI+1Zr2N2vkgw4NFjnvwdIOu5+orVmyEVZtKf7XUM34/3LFQD2lFUijzgT62xczarJ0PE3bDmVW16FnXRcnZQ2m5BOXtRVIn84E+tsXM2qydD3PT5efTP2N6AH378NHu5+mrRcrZg/L2IukTqbyy18WymFkMtfNBhgeL/Pl3dvLLd069iRw55p1f5KyRysJm4/c3bltpo8XQRLou8yP62MScn6928J2AvxRIUZ6+WtQHtKC8vUhKKNBHkVB+vqIn8vTVKg9o625cUqa8vUjXKdBHkVB+vqJn8vTVlqyGm18rj+4bUd5epJsyH+hHt0+ydMMzLFz3XZZueKa1wJlQfr5ieLDIme+e/rjkyDHniyMvpDfYQ1UqJ4Lx+2H9LKVzRDos04E+tuUPEszPV4Tl6VMzU7aeZvL2UE7n/LFG+CIdkulAH8vyBwnn5yvqTeBKzUzZeprJ21conSPSEZkO9LEsf5Bwfr4ibO2bisk0VuDUaipvX6Z0jkjiMh3o217+YGIEDu0JPhdTfr6i3to3UEqKpDp9U63ZVA6cTOf81VwFfJGYZTrQt7X8QWVJ4jAx5ucrhgeL/PXqD4RuwLt+bGfsPzMxraRyAA6/XQr4GuWLxCbTgb6t5Q/CliSGpteeb8bwYDF0V/WDU0d6Z1QPJ1M5q+6D/jOaf71G+SKxMA9ZPbGbhoaGfHx8vLudWD8r/Nyq+xLdWWnphmdCc/IDhX523PrxxH52oh6/sbQNYeitLILCOaW9bbWzlcgpzGybuw8Fncv0iL7lGvqJEULzy7PmJx5k6qWWem5UX63VdE61yih//Sz4m9/RSF8kgswG+rZq6L93M8GjTkssZVNteLDI2aeHL3+Q+klU9bSbzql2aI/y+SIRZDbQt1xDPzFSGjUG8o6lDG698v2h546584WHd/Dl0Rc70pdELFkNX9pXCvjtjPArqkf6Cvwip8hsoG+5hv7p28LPzZrfRo+a02hU78A3t/y8t4M9nBzhrz8Uzyi/ojbwK/hLjmU20LdcQx9WNw8dSdtUu/XK99edRAXwrS0/7900Tq24R/m1FPwlpzJbdfPl0Rf51pafn5JpL/T3hZdXToyUcvNhaZvCOaWRZ4eNbp/kiyMvcKzB+1QcKHDT5eena6OSOEyMwHdugCNvd+5n2ozShu+z5pdu7qrwkR5Qr+omk4G+8iC2OkdvwH+55Dz+Ynjx9BdUJkeF1c1jpWqRLv2DH90+yRce3tGwKLHuNfa6iZFSWq3eX1ydpDJPSZncBfqwOvTiQIEfrLts+gvuWFjnAWzZ+kMt9ycOXx59kW9u+XmktmfM7OMvP9nCvri9otFfX2miG4J0SK4C/ej2SW54eEfgOQNe2/CJUw8+fmPjPVBnzYcvvNRSf+LUTLCvOPv0fm698v3ZDfrQW4E/Kt0gpEm5CfRBeflqp4zoowaH/gJceVdq/sEN3vbktI3Em5Xp4H9KisdoaxZululGkjm5CPSNRruF/j4e/PDP+PArG6KP/FL4jyHo+UOSZhgc9x5/2KvgL7VS+G+7XW0HejNbBvwt0Ad83d031Jy38vkrgHeAP3T356O8NkizgX50+yTP/p+/4yunPcg59taJ407VQgbW1KK5XauyiWJ0+yR3bt7VG2vUp9hVM57j1prfmYqQ1aIlQ9xrYkTleM2xKG2iHovS5qCdxe4L/ycfvupz9S+gRr1AP32j0ukv7gPuBj4G7AW2mtmYu79c1Ww5sKj8cTFwD3BxxNe2bcd3N/HV/v/Fu+zUUW7r/1atdLdPqeHB4omRdSt5eykZO34pY4cvnXY86AZQ+49RN4LeZyGDv9pjUdpEPRalzdm8yQe2/Q+2QtPBPkzDQA9cBOx291cBzOwhYCVQHaxXAg966c+DLWY2YGZzgAURXtu2Pzr8Td41I65UhsHQZ3vmT7q/GF7M0G+dw/qxnRycai93LyVhN4CKsL8EgkZruiFIK2baUeY/fyd0MNAXgeri5b2URu2N2hQjvhYAM1sDrAE477zzInTrpLkz/l9T7UP1aN6ueoRfyuFPMHXkeJd7lV2NbgQVUW8IukFIkPf4L2L7b0UJ9GEbHkVpE+W1pYPum4BNUMrRR+jXCf9eeC+nT+1v5iXTDV1bWka3x1WCfnUeX48fuyPqDaFWM38xtJpP1o0k/V63c3lvTP+tKIF+L1C9mtc8YF/ENjMjvLZtpy+/jWOP/Ql9frT5F/foKL6R6lF+hYJ/b2j1BhFVlGcQ7RxLw3+rXpteuMkd9tPY86GbYgv0DatuzOw04F+AjwKTwFbg9919Z1WbTwDXU6q6uRi4y90vivLaIC3V0Ueti89oYE/K6PZJ5f8lM2pvcmm8cSVRdRO1vPIK4GuUSiQfcPe/NLO1AO5+b7m88u+AZZTKKz/j7uNhr23081KxlaCISA/JxYQpEZE8y+2esSIiokAvIpJ5CvQiIhmnQC8iknGpfBhrZgeAnwHnAvFND+s9eb5+XXt+5fn627n233L32UEnUhnoK8xsPOwpch7k+fp17fm8dsj39Sd17UrdiIhknAK9iEjGpT3Qb+p2B7osz9eva8+vPF9/Itee6hy9iIi0L+0jehERaZMCvYhIxqU20JvZMjPbZWa7zWxdt/uTNDP7qZm9aGY7zKyy8uc5ZvaUmf2k/PnsbvczDmb2gJm9bmYvVR0LvVYzu6X8e7DLzC7vTq/jE3L9681ssvz+7yiv+lo5l5nrN7P5Zvasmb1iZjvN7PPl45l//+tce/Lvvbun7oPSksb/Cvw2pc1LXgAu6Ha/Er7mnwLn1hz7KrCu/PU64I5u9zOma/1d4ELgpUbXClxQfv/fBSws/170dfsaErj+9cB/C2ibqesH5gAXlr8+i9J+FRfk4f2vc+2Jv/dpHdGf2JDc3Q8DlU3F82Yl8I3y198AhrvXlfi4+/eB2h1iwq51JfCQu//a3V8DdlP6/ehZIdcfJlPX7+773f358tdvAq9Q2ls68+9/nWsPE9u1pzXQh202nmUOPGlm28obpQP8prvvh9IvCfCervUueWHXmqffhevNbKKc2qmkLjJ7/Wa2ABgEfkjO3v+aa4eE3/u0BvrIm4pnyFJ3vxBYDlxnZr/b7Q6lRF5+F+4B/iPwQWA/8Nfl45m8fjM7E3gEuMHdf1WvacCxnr7+gGtP/L1Pa6CPsiF5prj7vvLn14HHKP2J9m9mNgeg/Pn17vUwcWHXmovfBXf/N3c/5u7Hgfs4+Sd65q7fzPopBbpvufuj5cO5eP+Drr0T731aA/1WYJGZLTSzmcDVwFiX+5QYMzvDzM6qfA18HHiJ0jV/utzs08A/dqeHHRF2rWPA1Wb2LjNbCCwCftSF/iWqEuTKPknp/YeMXX95f+n7gVfcfWPVqcy//2HX3pH3vttPous8ob6C0lPpfwW+1O3+JHytv03p6foLwM7K9QK/ATwN/KT8+Zxu9zWm6/0HSn+iHqE0arm23rUCXyr/HuwClne7/wld//8GXgQmyv/A52Tx+oFLKaUfJoAd5Y8r8vD+17n2xN97LYEgIpJxaU3diIhITBToRUQyToFeRCTjFOhFRDJOgV5EJOMU6EVEMk6BXkQk4/4/EG7ZSAIjdvQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(nospin['Bamb'],nospin['normalized'])\n",
    "plt.scatter(yesspin['Bamb'],yesspin['normalized'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
